close all
clear
%clc
addpath(genpath(pwd))
warning off
%% Section 1: User settings for the estimation
startYear     = 1961;       %Start year for estimation
endYear       = 2019;       %End year for estimation

lambdaSMM     = 1           %Weight in objective function to Shrinkage of based on Method of Moments (SMM)
SMMwithCSreg  = 1           %2 for incl. CS and risk-adj CS moments,
                            %1 for incl. CS moments,
                            %0 for not including CS moments
inclSurveys   = 0;          %1 for including surveys, else 0

appOrder      = 2;          % Approximation order of DSGE model
appOrderVar   = [appOrder;appOrder;appOrder]; % Used appOrder per variable: c_t, pi_t, evf_t
endoMeanGrid  = 0;          % 1 for endogenously recenter the grid given the projection approx.
                            % 0 for using the mean in Perturbation to re-center the grid 

muGrid        = 2          % mu in Smolyak grid
gridFactor_x  = 1.5         % grid for proj using +/- gridFactor_x*std1st approx
nx            = 6;          % #of states
ne            = 5;          % #of shocks

numCPUs       = feature('numcores'); % Number of CPUs

numOptimSteps = 0;          % Number of optimizations
optim         = 1;          % 1 Nelder-Mead.
                            % 2 Nelder-Mead with transformed parameter space
                            % 3 Own gradient based optimizer
                            % 4 CMAES
                            % 5 fmincon
cmaesPopSize  = 48;         % population size for CMAES
cmaesSigma    = 0.01;       % sigma in CMAES
cmaesSigmaMax = 0.05;       % sigmaMax in CMAES
MaxIter       = 50000;      % Max number of iterations for the optimizer
MaxEvals      = 10000;      % Max number of function evaluations for the optimizer
TolFun        = 1D-3;       % Function tolerance at the objective function
TolX          = 1D-3;       % Function tolerance for the coefficients
epsValue      = 1e-5;       % stepsize for computing gradient in own optimization routine
seOn          = 1;          % 1 for computing standard errors, else 0
ScalingMacro  = 0.2;        %Setting measurement errors for macro variables
ScalingYield  = 0.1;        %Setting measurement errors for yields

if numCPUs > 8
    % We use Grendel-S
    rmpath([pwd,'/Mex_myPC'])
    MEXon     = 1;          %1 For using MEX-files, else 0
else
    % We use My PC
    rmpath([pwd,'\Mex_GrendelS'])
    MEXon     = 1;
end

% Default setting
tau           = 400*0;      %Moments computed based on T*tau simulations, else closed-form solution used when tau = 0.
pruningOn     = 1;          %1 for using the pruned approximation, 0 for the unpruned approximation
nyMom         = 10;         %number of the first nyMom observables used for the SMM part
CSregress_m   = 4;          %2 for biannual implementation of CS regressions, 4 for annual etc.

%% Section 2: Data and Model parameters
% Setting the calibrated parameters as in the sample ending in 2007Q4
structData = loadData_2020(startYear,endYear);
meanData   = nanmean(structData.data,1)';
stdData    = nanstd(structData.data,0,1)';
setCalibrateParams

CSdata     = CampbellSchillerRegression(structData.yieldsAll(1:end,:),CSregress_m);
%[CSdata.maturities;CSdata.CSbetta(2,:); CSdata.CSBetta_CI95Boot]
%[CSdata.maturities;CSdata.CSbetta(2,:);(CSdata.CSbetta(2,:)-1.96*CSdata.CSbetta_se(2,:));(CSdata.CSbetta(2,:)+1.96*CSdata.CSbetta_se(2,:))]
PredData    = predictiveRegression(structData.yieldsAll(1:end,:),CSregress_m);


% Compute empirical moments and weigthing matrix
dataSMM = momentsSMM(structData.data,structData.yieldsAll,nyMom,CSregress_m,SMMwithCSreg);
qLag = 3;
W    = getOptimalWeighting(qLag,dataSMM);

% The variables which we do not select for estimation are calibrated to the assigned values
%calibrateParams        = rmfield(calibrateParams,'U0d'); %use this to have U0d being active
%calibrateParams        = rmfield(calibrateParams,'U0'); %use this to have U0 being active
%calibrateParams.CRRA    = 10;
calibrateParams.RHOz     = 0;

% Starting values
params.ALFA    = -20;
params.BETTA   = 0.995;
params.B       = 0.1;
params.PHI     = 0.075;
params.CHI     = 2.8;
params.XI      = 0.75;
params.KAPAw   = 0;
params.ETA     = 6;
params.PHIpai_1= 2;
params.PHIc_1  = 0.0;
params.RHOz    = 0.0;
params.RHOd    = 0.99;
params.RHOn    = 0.95;
params.RHOp    = 0.99;
params.RHOa    = 0.99;
params.PAIss   = 1.0;
params.SIGMAmuz= 0.0004;
params.SIGMAd  = 0.03;
params.SIGMAn  = 0.03;
params.SIGMAp  = 0.001;
params.SIGMAa  = 0.01;

% With select = [2:10]*4/CSregress_m in SMM
if lambdaSMM  == 1 && SMMwithCSreg  == 1
    calibrateParams.ScalingYield = 0.083;   %corresponds to 25bp for yields and surveys
    
    paramsValues = [-16.48753574  0.99336222  0.86609872  0.20319048  1.94827519  ...
        0.77195767  0.86399910  8.15706975  2.94492588  2.33383626  ...
        0.76954953  0.97798508  0.96661087  0.99653933  0.99030401  ...
        1.01834696  0.00297334  0.03329981  0.03502314  0.00082606  ...
        0.00392939 ]';
    % Q =   -11.5152, done 

    if inclSurveys   == 0  
        paramsValues = [-16.488142627127    0.993361796882    0.866151653335    0.208673821947    1.948284600473 ...
            0.771988917925    0.866705114869    8.208260170706    2.981735735538    2.343907644462...
            0.770458142825    0.978012102962    0.966639957868    0.996530002460    0.990317211615...
            1.018413492230    0.002973370015    0.033300207827    0.035043870760    0.000826204153...
            0.003954203540 ]';
        % Q =    -1.2318, done
    end


end

selectParams = fieldnames(params);
params       = values2struct(paramsValues,selectParams);

%% Section 3: Setting setupFilter and setupProj
% The lower and upper bounds for the parameters
[lowerBounds,upperBounds,Insigma] = setLowerUpperBounds;

% The information needed for projection
Smolyak_elem_iso                  = Smolyak_Elem_Isotrop(nx,muGrid);
setupProj.Smol_grid_iso           = Smolyak_Grid(nx,muGrid,Smolyak_elem_iso);
[n_nodes,epsi_nodes,weight_nodes] = Monomials_2(ne,eye(ne));

% Constructing the struct setupFilter
constructSetupFilter
setupFilter.zlbON          = 1;
%% Section 3: Estimation
% tic
 [Q0,errorMes,model0,~,outFilter0] = objectFuncQML(struc2values(params,setupFilter.selectParams),setupFilter);
% Q0
% toc

if numOptimSteps > 0
    for i=1:numOptimSteps
        % We switch between the optimizer in optim and the Nelder Mead
        if mod(i,2) == 1
            setupFilter.optim = optim;
        else
            setupFilter.optim = 6;
        end
        outEstim = runOptimization(setupFilter,params);
        params = outEstim.paramsOpt;
        printToScreen(struct2array(params))
    end
else
    setupFilter.optim = 0;
    outEstim = runOptimization(setupFilter,params);
end

%% Model diagnostics
outEstim.outFilter.Q
outDiagnostics = plotOutputFilter(outEstim);
plotTermPremia(outEstim)
plotCSloadings(CSdata,outEstim);
reportFit
CRRA = getCRRA(outEstim.model.params)
outEstim.newsDecomSim = inflVarianceRatio(outEstim,outEstim.outFilter.outSim.xAll_cu);
outEstim.newsDecom    = inflVarianceRatio(outEstim,outEstim.outFilter.xhat);
outEstim.newsDecomSim.summaryTable
%outEstim.newsDecom.summaryTable
outEstim.termPremSim.stdTP(end)
%% Additional ex post estimation analysis
outEstim.UnconditionalMoments = unconMoments(outEstim.outFilter.ySimAsInData,outEstim.setupFilter);
outEstim.reducedForm          = doReducedFormAnalysis(outEstim);

% Computing GIRFs
close all
IRFlength = 20;
shockSizeScalar = 1;
xsGIRF             = [mean(outEstim.outFilter.xhat(7,:)); zeros(5,1)];
outEstim.GIRFs     = getImpulseResponses(outEstim,IRFlength,shockSizeScalar,[],xsGIRF);

plotGIRFsChart(outEstim.GIRFs,2)

xfGIRFxhat         = outEstim.outFilter.xhat(1:6,:);
xsGIRFxhat         = [outEstim.outFilter.xhat(7,:);zeros(5,setupFilter.T)];
outEstim.GIRFsxhat = getImpulseResponses(outEstim,IRFlength,shockSizeScalar,xfGIRFxhat,xsGIRFxhat);
xfGIRFsim          = outEstim.outFilter.outSim.xAll_cu(1:6,1:1000);
xsGIRFsim          = [outEstim.outFilter.outSim.xAll_cu(7,1:1000);zeros(5,1000)];
outEstim.GIRFsSim  = getImpulseResponses(outEstim,IRFlength,shockSizeScalar,xfGIRFsim,xsGIRFsim);

plotGIRFsChart(outEstim.GIRFsSim,2)

testOfSpanningON = 1;
if testOfSpanningON == 1
    numPCAs             = 5;
    outEstim.macroSpanData_5 = doMacroSpanRegressions(structData.yieldsAll(1:end,:),structData.data(1:end,:),setupFilter,numPCAs);
    outEstim.macroSpanData_6 = doMacroSpanRegressions(structData.yieldsAll(1:end,:),structData.data(1:end,:),setupFilter,numPCAs+1);
    % On Simulated data with Measurement errors
    yieldsAllSimWithNoise        = outEstim.outFilter.yieldsAllSimWithNoise;
    ySimAsInDataWithNoise        = outEstim.outFilter.ySimAsInDataWithNoise;
    outEstim.macroSpanSim_5      = doMacroSpanRegressions(yieldsAllSimWithNoise,ySimAsInDataWithNoise,setupFilter,numPCAs);
    outEstim.macroSpanSim_6      = doMacroSpanRegressions(yieldsAllSimWithNoise,ySimAsInDataWithNoise,setupFilter,numPCAs+1);
    
    % On Simulated data withOUT Measurement errors
    yieldsAllSim = outEstim.outFilter.yieldsAllSim;
    ySimAsInData = outEstim.outFilter.ySimAsInData;
    outEstim.macroSpanSimNoNoise_5 = doMacroSpanRegressions(yieldsAllSim,ySimAsInData,setupFilter,numPCAs);
    outEstim.macroSpanSimNoNoise_6 = doMacroSpanRegressions(yieldsAllSim,ySimAsInData,setupFilter,numPCAs+1);
    
    
    % a linear model - without measurement noise
    modelLin         = outEstim.model;
    modelLin.gxx     = modelLin.gxx*0;
    modelLin.gxxx    = modelLin.gxxx*0;
    modelLin.Gxx     = modelLin.Gxx*0;
    modelLin.Gxxx    = modelLin.Gxxx*0;
    modelLin.hxx     = modelLin.hxx*0;
    modelLin.hxxx    = modelLin.hxxx*0;
    modelLin.Hxx     = modelLin.Hxx*0;
    modelLin.Hxxx    = modelLin.Hxxx*0;
    modelLin.byxx    = modelLin.byxx*0;
    modelLin.Byxx    = modelLin.Byxx*0;
    modelLin.byxxx   = modelLin.byxxx*0;
    modelLin.Byxxx   = modelLin.Byxxx*0;
    if inclSurveys  == 1
        modelLin.Rxx     = modelLin.Rxx*0;
        modelLin.Rxxx    = modelLin.Rxxx*0;
    end
    setupFilterLin   = outEstim.setupFilter;
    setupFilterLin.pruningOn = 0;
    
    % Moments computed by simulation
    setupFilterLin.numSim = 1D5;
    outSimLin      = simProjection(modelLin,setupFilterLin.shocks(:,1:setupFilterLin.numSim));
    if setupFilter.pruningOn == 1
        yieldsSimLin  = getYields(modelLin,outSimLin.xAll_cu);
    else
        yieldsSimLin  = getYields(modelLin,outSimLin.x_cu);
    end
    pos_l             = find(strcmp({'$l_t$'},modelLin.labely));
    pos_w             = find(strcmp({'$w_t$'},modelLin.labely));
    pos_c_ba1         = find(strcmp({'$c_{t-1}$'},modelLin.labelx));
    pos_myz           = find(strcmp({'$\mu _{z,t}$'},modelLin.labelx));
    pos_c             = find(strcmp({'$c_t$'},modelLin.labely));
    pos_pai           = find(strcmp({'$\pi_t$'},modelLin.labely));
    ySimAsInDataLin   = zeros(setupFilterLin.numMacro+length(modelLin.maturities),setupFilterLin.numSim);
    ySimAsInData(1,:) = outSimLin.y_cu(pos_l,:)-modelLin.gSS(pos_l,1);
    ySimAsInData(2,:) = outSimLin.y_cu(pos_w,:)-modelLin.gSS(pos_w,1);
    ySimAsInData(3,:) = 4*(outSimLin.y_cu(pos_c,:)-(outSimLin.x_cu(pos_c_ba1,:)+modelLin.gSS(pos_c_ba1,1))+outSimLin.x_cu(pos_myz,:)+log(modelLin.params.MUZss));
    ySimAsInData(4,:) = 4*outSimLin.y_cu(pos_pai,:);
    ySimAsInData(setupFilter.numMacro+1:setupFilter.numMacro+length(modelLin.maturities),:) ...
                         = 4*yieldsSimLin(modelLin.maturities,:);
    yieldsAllLin        = 4*yieldsSimLin;

    outEstim.macroSpanSimNoNoiseLin_5 = doMacroSpanRegressions(yieldsAllLin,ySimAsInData,setupFilter,numPCAs);
    outEstim.macroSpanSimNoNoiseLin_6 = doMacroSpanRegressions(yieldsAllLin,ySimAsInData,setupFilter,numPCAs+1);
end

% Deleting some simulated series
outEstim.outFilter = rmfield(outEstim.outFilter,'yieldsAllSim');
outEstim.outFilter = rmfield(outEstim.outFilter,'yieldsAllSimWithNoise');
outEstim.outFilter = rmfield(outEstim.outFilter,'ySimAsInData');
outEstim.outFilter = rmfield(outEstim.outFilter,'ySimAsInDataWithNoise');

% Standard naming of output file
filename = ['Model_endYear',num2str(endYear),'_lambdaSMM',num2str(lambdaSMM),'_SMMwithCSreg',num2str(SMMwithCSreg),'_inclSurveys',num2str(inclSurveys)];
if seOn == 1
    filename = [filename,'_withSE'];
end
save(filename,'outEstim','outDiagnostics','structData','CSdata')
%[struct2array(outEstim.paramsOpt)' struct2array(outEstim.paramsSE_eps6.qLag3.paramsSE)' struct2array(outEstim.paramsSE_eps7.qLag3.paramsSE)' ]